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Abstract 

The purpose of this paper is to describe the application of the notion of viscosity solutions to solve the 
Hamilton-Jacobi-Bellman (HJB) equation associated with an important class of optimal control problems 
for quantum spin systems. The HJB equation that arises in the systems of interest is a first order nonlinear 
partial differential equation defined on a Lie group. We employ recent extensions of the theory of viscosity 
solutions from Euclidean space to Riemannian manifolds to interpret possibly non-differentiable solutions to 
this equation. Results from differential topology on the triangulation of manifolds are then used to develop 
a finite difference approximation method, which is shown to converge using viscosity solution methods. An 
example is provided to illustrate the method. 

Keywords: optimal control, quantum spin systems, dynamic programming, numerical method, viscosity 
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1. Introduction 

Recently, there has been considerable attention 
directed at the problem of obtaining time optimal 
trajectories for open loop control of quantum spin 
systems [1, 2, 3, 4]. These problems arise from 
applications which include NMR spectroscopy (to 
produce a time optimal trajectory), and the opti- 
mal construction of quantum circuits [5, 6] (to mini- 
mize the number of logic gates required to construct 
a desired unitary transformation). These spin sys- 
tems have the mathematical structure of a bilinear 
right invariant system on the special unitary group. 
Owing to the importance of the applications, there 
have been various approaches to solving these prob- 
lems which utilize Lie theoretic arguments [2, 7], 
calculus of variations [4, 8] and Dynamic Program- 
ming [9, 10]. 

In the dynamic programming approach, under 
appropriate regularity assumptions, the optimal 
cost function (value function) is the solution to a 
Hamilton- Jacobi-Bcllmann (HJB) equation [11, 12, 
13]. For many problems of interest this value func- 
tion can be demonstrated to be non-diffcrcntiable. 
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Hence there is the need for a more general notion 
of a solution to such PDEs. A popular and suc- 
cessful concept of such a weak solution of nonlinear 
PDEs is the well studied theory of viscosity solu- 
tions [14, 15] on Euclidean spaces. Because the 
quantum spin problem leads to a HJB equation 
defined on a Lie group, we use extensions of the 
viscosity solution theory to Riemannian manifolds 
[16, 17, 18] in order to interpret the solutions of 
this equation. For a detailed introduction to this 
topic we refer the reader to [14, 19] and references 
contained therein. 

In this article we build up the components re- 
quired for a rigorous application of viscosity so- 
lution theory on manifolds for quantum systems. 
This commences with an explanation of a dis- 
cretization method based on the triangulation of 
manifolds [20] to solve the HJB equation for the 
optimal spin control problem. We then use viscos- 
ity solution concepts to prove the convergence of 
the solution obtained by this triangulation-based 
discretization scheme to the solution of the original 
HJB equation. 

The structure of this article is as follows. We 
begin by describing the quantum spin control prob- 
lem in Sec 2. This is followed in Sec 3 by a study of 
the regularity properties of the value function which 
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play an important role in the solution of the asso- 
ciated HJB equation. After motivating the need 
for more generalized solutions of the HJB equation 
using an example system with a non-diffcrentiable 
value function, we explain the use of the notion of 
viscosity solutions on Lie groups in Sec 4. Results 
pertinent to the existence and uniqueness of such 
solutions are recalled from relevant literature and 
modified to the framework of the problems intro- 
duced. In order to solve these optimal control prob- 
lems numerically we make use of the notion of trian- 
gulation of the group on which the system evolves. 
This concept and the proofs of convergence of the 
approximations to the actual solution using viscos- 
ity solution notions are introduced and developed 
in Sec 5. In Sec 6 the ideas developed are then ap- 
plied to solve an example control problem on SU (2) 
for which sample optimal trajectories and the value 
functions from the simulations arc obtained. We 
conclude with comments and possible extensions in 
Sec 7. 



by U(t;v,Uo). We denote by T a compact set in G 
with smooth closed boundary dT. This set is the 
target set that we wish the system to reach. 

Now, the following properties are satisfied at ev- 
ery point x e G: 

• the system dynamics /(•) is driftless, 

• if X := f(x,u) can be generated by a certain 
control u e V, then —X can also be generated 
using another element of the control set (in this 
case it would be — u), 

• the dimension of the vector space at x, gener- 
ated by the set of vector fields after all possible 
bracketing operations, has the dimension equal 
to that of the group. 

Hence we have from [21, Prop. 3.15] that the time 
to get from any point on the group, to the iden- 
tity element is bounded. Thus the entire group is 
reachable from the identity and hence the problems 
dealt with in the next section are well defined. 



2. Problem Description 

2.1. System Description 

In this section we introduce a mathematical 
model arising in applications of the open loop con- 
trol of quantum systems. Given a compact con- 
nected matrix Lie group G = SU(2 n ) with an as- 
sociated Lie algebra g — su(2") and smooth, right 
invariant vector fields X\, ....X m in q, let the evo- 
lution of the system be given by 



dU 
~dt 

U,U e G. 



= {Y / v k (t)Xk}U, U(0) = U , 

(1) 



fe=i 



Here X\ , . . . , X m are termed control vector fields. 
The Vi are elements of the control signal v which be- 
longs to the class of piecewise continuous functions 

V with their range belonging to a compact subset 

V of the real m-dimensional Euclidian space K m , 
containing the origin. Without loss of generality 
we may consider V to be the unit hypercube in R m 
around the origin. We assume that the Lie algebra 
generated by the set {X\, X m }, using repeated 
Lie bracketing operations of all orders, is g. We de- 
note the right hand side of the Equation (1) above 
by f(U,v). Given a control signal v and an initial 
point Uq, the solution to Eq (1) at time t is denoted 



2.2. Problem Formulation 

A large class of problems in the control of quan- 
tum spin systems and quantum circuit synthesis 
can be recast in terms of an optimal control prob- 
lem with the following value function (optimal cost 
function) : 

!tu (v) 
J £(U(s;v,U ),v{s))e- Xs ds 

where t : G x V — > R is continuous and A > is 
a real valued discount factor. Here t Uo (v) denotes 
the time to reach the set dT starting from Uq using 
control v. Note that under the assumption that the 
cost £(.) has a lower bound which is positive, we 
can then assume without loss of generality that the 
minimum of I over G x V is 1. We now proceed to 
study some properties of the value function defined 
above. 

3. Regularity of the Value Function 

We begin by introducing certain quantities which 
will be used in the study of the regularity properties 
of the value function. Define the minimum time 
function 



T(U) :- 



M{tu(v)}, 



u e G, 



(3) 
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which is the infimuni of the time taken to reach 
the target set from a starting point U. The set of 
points, termed the reachable set, from which the de- 
sired target set may be reached in time r is defined 

as 

R(t) :={U e G | T(U) <t}. 

A system is said to be small time controllable on a 
set T (denoted by STCT) if 

TCint(R(t)), Vt>0. (4) 

Note that the assumption of small time controlla- 
bility to the target set is required to obtain some of 
the results in this section and holds for the exam- 
ples studied in this paper (and may be verified for 
any system under consideration). 

We now introduce some results on the regularity 
of the value function whose proofs proceed along 
the lines of the arguments used in [14] with suitable 
modifications due to the Lie group setting. 

Lemma 1 ( Proposition 1.2 (§ 4) [14]). If the 

system is small time controllable on the set T then 
the value function is continuous in some open set 
containing the boundary dT of the set. 

Lemma 2. Given a system evolving on a con- 
nected, compact Lie group, with dynamics ( Eq (1) ) 
such that f , I, S satisfy the following conditions: 

1. / : G x V — » TG is continuous. 

2. <! : G x V ^ R is continuous. 

3. S is continuous on some open set containing 
dT 

Then S is bounded and continuous on G. 

Proof. This result consists of two cases of the 
value function. 

• A = 0: As indicated in Sec. 2.1, the time to get 
from any point on the group, to the identity el- 
ement is bounded. This along with the bounds 
on £, imply that the value function is bounded. 

• A > 0: In this case the boundedness follows 
directly from the bounds on I and the expo- 
nential decay factor A > 0. 

In both cases, the continuity proofs proceed along 
the lines of [14, Prop 3.3 (§ 4) ] with suitable mod- 
ifications for the Lie group setting. 

Example 3 (Property of Value Function). 



Consider the following system defined on the special 
unitary group G = SU(2): 

U=[I X + vl x ] U, (5) 
U(0) = U , 

where v : [0 , oo) — > (—00,00) is a piecewise con- 
tinuous control signal. Here I x and I z are given 
by: 

2 V 1 J ' 2 ^ -1 ) • 

(6) 

Note that in Eq (6), we follow the mathematicians' 
convention in which I x and I z are skew-Hermitian 
matrices which belong to the Lie algebra su(2) of 
the group SU(2). The cost function to reach the 
target set for this problem (which in this case is 
the identity element J) is given by 

tu (v) 

S{U )=ini I cxp{-A s}ds. (7) 
vev J 


We demonstrate below that this normalized mini- 
mum time function S : SU{2) — > R is not differ- 
entiate at the points ±1. 

Any point P in SU{2) can be represented as: 
k\ exp(a I z ) &2 where k\ , are elements of the Lie 
subgroup generated by exp(I x ) and a £ [0 , n). Us- 
ing [1], the cost function for any point expressed in 
this form is given by 

1 — exp{— Aa} 
A^ ' 

For this compact, connected Lie group the expo- 
nential mapping (denoted by (p) is a diffcomorphism 
from an open set around the origin in su(2) to an 
open set around / in SU(2). Let the I z axis in 
su(2) be denoted by e\. If the function S is diffcr- 
cntiable at the identity element then the function 
S := S o : su(2) — > R must be differentiate at 
the origin in su(2). Hence there must must exist a 
linear function 77 : M 3 — > M s.t 

hm l^ + e)-S( X )-r ]{ e)\ at a; = 0. 

IMKo [ ||e|| J 

(8) 

Now, consider a line through the origin in su(2) 
along the e\ axis. Let e be either +8e\ or — 8e\ 
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(with S > 0). The value of S has the following 
properties at the origin 

5(0) = 

S( + Se 1 ) = S(-Se 1 )= 1 -^^l. (9) 

In addition, at the origin the function 77(e) takes the 
form SL where 

L := ei • DS(0), 

where DS denotes the differential of the value func- 
tion. Hence using this expression for 77(e) and 
Eq (9) in Eq (8), it follows that for S to be dif- 
fcrcntiable at the origin 



lim 



j S(±5. 



SL 



0. 



(10) 



Thus from Eq (9) and the Taylor expansion for 
cxp{A<5} it follows that the two equations below 
must simultaneously hold in order to ensure dif- 
ferentiability 



lim (s-SL + o^y 

8^0 { s 

<5->0 5 



(11) 
(12) 



The first of these requires L = +1 and the second 
requires L = — 1, leading to a contradiction. Hence 
S is not differcntiable at the origin, 

and therefore the function 5* is not differentiablc 
at /. Similar arguments hold for the element —I of 
SU{2). □ 

In view of this example, care is needed when in- 
terpreting the notion of solutions to the HJB equa- 
tion since the value functions may not be diffcr- 
entiable. This motivates our study of the use of 
non-differentiable weak solutions (more specifically 
viscosity solutions) to the HJB equation [14]. 

4. Viscosity Solution Theory 

4-1. Introduction 

The example in the previous section indicates 
the need for notions of weak solutions to the HJB 
equation associated with the dynamic programming 
problem. In this section we recall certain recent ex- 
tensions of viscosity solution theory to Riemannian 
manifolds [16, 17, 18] and present them in a form 



suitable to our study of systems evolving on the 
special unitary group G. We then indicate certain 
results about the existence and uniqueness of such 
solutions. 

Apart from being useful as a weak solution to 
the HJB which we require, the viscosity solution 
concept will also be used to prove convergence (in 
a specific sense) of the numerical approximation 
schemes to the original problem. The definitions 
and results in this article are more generally appli- 
cable to Riemannian manifolds, but for our current 
problem formulation the Lie group setting is suf- 
ficiently general. Note that we use the notation 
C k (Q) for the class of functions whose fc-th order 
derivatives are continuous on a set fi. In what fol- 
lows Dcf> : TG — > M. denotes the differential of a 
function cf> : G -> R. 

We start by defining the notion of viscosity solu- 
tions from [14, Chapter 2] 

Definition 4 (Continuous Viscosity Solution). 

Given an open domain ft in a Lie group G, a 
function S € C(f2) is a viscosity sub(super) 
solution of the following PDE in D, 

F(U,S,DS) =0, U e n 

ifV<p e C^H) 

F(U,S,D(f>) < (>)0 

at every point U <E where S — (j) has a relative 
maxima (resp. minima). A function is a viscos- 
ity solution iff it is both a super and sub viscosity 
solution. 

4-2. On the Existence and Uniqueness of Viscosity 
Solutions 

We proceed to develop certain results regarding 
the viscosity solution to the HJB equation arising 
from the associated control problem. 

Lemma 5. The value function S given by Equa- 
tion (2) is a continuous viscosity solution of the 
HJB equation 



F{U, S, DS) :=\S + H(U, DS) = 0, 
U e G\T, A > 



(13) 



where 



H{U,DS) := 



sup {-DS(U) ■ f(U,v)-t(U,v)}. 
vev 
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PROOF. This follows along the lines of [14, Prop 
3.11 (§ 4)], with slight modifications. 

Now, in order to prove uniqueness of the viscosity 
solution, we recall the following lemma. 

Lemma 6 (Corollary 4.10 of [17]). Let G be a 

compact Lie group and ft be a bounded open subset 
of G. Suppose that F : T*G x R -» R is uniformly 
continuous. Given that Vr > s 

3j > s.t 7(r -s)< F(x, r, Q - F(x, s, (), 

(14) 

where x G G. There is at most one viscosity solu- 
tion of the HJB equation 

F(x,S,DS(x)) = 0. 

We arrive now at the main uniqueness results for 
the discounted value function. 

Corollary 7. Assume that the value function S 
(with A > 0) is continuous. It is the unique vis- 
cosity solution to the HJB equation (13) on G \ T . 

Proof. From Lemma 5, S is a viscosity solution 
to the HJB equation (13). In addition, this HJB 
equation can be shown to satisfy the conditions in 
Lemma 6. Hence the statement of this corollary 
follows. 

We now deal with the case of the undiscounted 
cost, starting specifically with the minimum time 
function, since it is of special interest in applica- 
tions. 

Theorem 1. Assume that the minimum time 
function S ( with t = 1 and A = ) is continuous. 
It is the unique viscosity solution to the HJB equa- 
tion (13) on G\T. 

Proof. In this case we have that S is a viscosity 
solution of the HJB equation 

sup {-DS(x) ■ f(x, v)-l} = 0, x e G. (15) 
vev 

Applying a particular diffeomorphism called the 
Kruskov transform to S we obtain a new value func- 
tion R as follows 

R(x) := l-e- s( - x \ (16) 

From the transform theorem for viscosity solutions 
[14, Prop 2.5, Ch.2] it follows that S is a viscosity 



solution of Eq (15) iff the function R is a viscosity 
solution of the equation 

sup {-DR{x) ■ f(x, v)-l + R(x)} = 0, (17) 
t>ev 

<=> R(x) + sup {-DR(x) ■ f(x,v) -1} = 0. (18) 

As R is continuous (due to the continuity assump- 
tion on S) , from Lemma 5 we have that R is a vis- 
cosity solution to the HJB equation (18). Moreover, 
this new HJB equation, denoted by F(x,DR,R) = 
(say) , has the F expression satisfying the require- 
ments of Theorem 6. Hence R is a unique viscos- 
ity solution to Eq (18) with boundary condition 
R(x) = 0,Va; e T. As the Kruskov transform 
is a diffeomorphism, it follows that the minimum 
time function is the unique viscosity solution to the 
HJB Eq (15). 

The corresponding result for an undiscounted 
cost function with a generalized running cost £(x, v) 
is now described. 

Theorem 8 (Existence and Uniqueness). 

Assume that the value function S (with A = 0) is 
continuous in G \ T ■ In addition assume that 

1. I : G x V -^R is uniformly Lipschitz continu- 
ous and bounded. 

2. / : G x V — > TG is uniformly Lipschitz con- 
tinuous, bounded. 

then S is the unique viscosity solution of 

H(x,DS) = 0, mG\T, (19) 

and R(x) = 1 — e -5 ^ is the unique viscosity solu- 
tion of 

R(x) + sup {-DR(x) ■ f(x, v) - £(x, v)} , m G \ T. 

ti£V 

(20) 

Proof. By a generalization of [14, Ch.4, Prop. 
3.12] to the Lie group setting, R is a viscosity solu- 
tion of Eq. (20) and, using the transformation the- 
orem for viscosity solutions [14, Prop 2.5, Ch.2], S 
is a viscosity solution of Eq. (19). The uniqueness 
follows as the forms of these Hamiltonians satisfy 
the conditions in Lemma 6. 

5. Approximating the Viscosity solution on 
Lie groups via Triangulations 

In this section we describe the main focus of this 
work on numerical methods based on discretizing 
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the HJB equation. This requires a discretization of 
the Lie group on which the system evolves. The 
intuitive idea is to obtain a grid on the Lie group. 
This is a generalization of the idea of triangulating a 
2-dimcnsional surface. Such triangulations are used 
in numerical approximation procedures to obtain a 
solution to a discretization of the problem [22] . In 
the latter part of this section we apply the notion 
of viscosity solutions to the HJB equation, to prove 
the validity of such numerical approximations. 

5.1. Setup for Discretization 

In order to describe the discretization we recall 
some definitions from [20, Section 7.1, Ch.II ]: 

Definition 9 (Simplex). If vo, Vi ■■■v m are in- 
dependent points of E™, the simplex a which they 
span is the set of points x such that x = S&j^j. 
Note that bi > and S6j = 1. The numbers bi are 
termed the barycentric coordinates of x and, due to 
the properties listed here, may be used as transition 
probabilities in the algorithms described below. 

Definition 10 (Face of a Simplex). A face of 
the simplex a is the simplex spanned by a subset 
of the vertices of a. 

Definition 11 (Simplical Complex). A simpli- 
cal complex K is a collection of simplices in E™ such 
that 

1. Every face of the simplex ofK is in K. 

2. The intersection of two simplices ofK is a face 
of each of them. 

3. Each point of\K\ has a neighborhood intersect- 
ing only finitely many simplices ofK where |K| 
denotes the union of simplices ofK. 

Definition 12 (Star). If x is a point of |K|, the 
star of x in K is the union of the interiors of all 
the simplices a such that x lies in a. It is denoted 
by St(x,K). 

Let u! be a C r map. Given a point b in a we 
define the map dujb : o — > E" as 

dw b (x) := Du(b) ■ (x- b). (21) 

We now recall the definition of a triangulation [20, 
Section 8.3, Ch.II ] , which is the main concept 
required for the discretization of the group. 



g (algebra) 




Figure 1: Mapping (with overlaps) from the algebra to the 
group. 

Definition 13 (Triangulation). A C r map u> : 

K — » G is said to be an immersion if duj x : 
St(x, K) — > E" is injective for each x. Such an im- 
mersion which is a homeomorphism onto, is called 
a C r triangulation ofG. 

From [20, Theorem 10.6] it follows that every dif- 
ferential manifold has a triangulation. Additionally, 
since our problem framework is on a compact con- 
nected Lie group there is a set of natural imbed- 
dings from complexes in the Lie algebra q to the 
group via the exponential mapping. Note that com- 
plexes in the algebra may be easily obtained in this 
case, for instance via triangulations on the algebra 
(since the algebra in our case is isomorphic to some 
Euclidean space). Note that there may be overlaps 
between the images (on the Lie group) of these com- 
plexes. For instance we may have C r imbeddings of 
two complexes into G whose images overlap. It can 
be shown that by altering these imbeddings slightly 
their images can be made to fit together 'nicely' 
i.e intersect in a subcomplex after suitable subdivi- 
sions. The intuition behind this is conveyed in Fig 
1. This is an especially important concept since 
mappings such as the exponential map are not in- 
jective and carry multiple points in the algebra to 
the same point on the group. 

The following notations are used in this section: 

• L° := G \ T. 

• L := closure of L°. 
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Sh'- simplex on G with maximum distance be- 
tween points being ft. 



• L 



• dL°h- subset of s.t 

sup inf ||x - y\\ -> as h -> 0. (22) 

Thus <9L£ is the approximation to the bound- 
ary of L°. 

We recall that the system is denoted by 

x = f(x, v) , x e G 

where v is a control signal chosen from a compact 
topological space of controls V. 

Now consider the continuous viscosity solution V 
to the HJB equation on the set L. The discretized 
version of the solution is denoted by V h . In the case 
of the exponentially discounted (with discounting 
factor A), infinite time horizon control problem [14, 
22] the HJB equation is 

XV + sup{-DV(x) • f(x, v) - £(x, v)} = 0, (23) 

with boundary condition 
V(x) = 0, x e dT. 

In order to discretize the HJB equation we define 
the following terms 

£^):= sup || /(z,«) ||, (24) 

and At h (x) := J- y (25) 

and use the following discretization for DV ■ f 
DV-f(x,v) = 



1 

ft 



p(x,y\v)V(y)-V(x)yf(x,v). (26) 



Here T, h (x) is the set of vertices of the simplex 
which contains the point arrived at by flowing for 
time At h (x) along f(x,v) using a constant control 
signal v starting from the point x. This is depicted 
in Fig. 2. For each value of v the terms p(x,y\v) 
must satisfy 

p( x >y\ v ) = 1 , 

y eY, h (x) 



, U{Af(x):v,i 



Figure 2: Discretization using Barycentric coordinates. 



and can therefore be interpreted as transition prob- 
abilities. Hence they can be obtained naturally 
from the Barycentric co-ordinates on the complex 
as mentioned previously. Using this discretization 
in the HJB equation (23) and rearranging we obtain 



\V h {x) 



ievi h i ^ p( x >y\ v ) v (y) - v ( x ) • f( x > v ) 




(27) 



Dividing throughout by B(x) and rearranging, 
we get 



\V(x)h 



B(x) 



inf 



1^.1 [^2 p(x,y\v)V(y)-V(x) ■ f(x,v 




(28) 



which implies 



V(x) 1 + 



Aft 

Bjx) 

E ve s" {p( x > y\ v ) v (y)f( x > v )} 



= inf J 

+ V ^>{ i B(x) J + B(x) J ' 
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Therefore, 

V(x) (l + \At h (x)) (29) 
Eyes* {p(x,y\v)V(y)f(x,v)} 



inf ( 
V{x)(\ 



B{x) 



f(x,v) 
B(x) 



) +e(x,v)At h (x)Y 



We define a new set of transition probabilities 

fix v) 

p(x,y\v) =p(x,y\v) ' , Vy ^ x (30) 



p(:r,2;|i>) = f 



B(x) 
B{x) ■ 



(31) 



Hence the discretization of the HJB equation can 
be shown to satisfy 

V h (x) = 



t>ev 



™l\ ( E P(x,y\v)V h (y) 



y£~t h (a 



+^^W)X IT ^}, (32) 

with boundary condition 

V h (x) = 0, a: e dL° h , 

where 

S k W:=E l W[JW. (33) 

We rewrite Eq (32) in the more general form 

V h {x) = F h {V h ){x). (34) 

Hence it can be seen that V is the fixed point of 
an operator F . We assume that there exists a pa- 
rameter -fh e (0, 1] which is continuous with respect 
to h with the property that lim 7/, = 1. It can be 

checked that F h (-) has the following properties: 

1. For all constant valued functions c 

F h (V h +c)(x) = F h (V h )(x) +1 ^c, x e K- 

2. ForallweC^L) 

Um F h (u)(y)-u,(y) 



y—*x 
hi.0 

lim 



> <=> H{x,Duj 1 uj) < 



F h (LQ)(y)-Lo(y) 



<0^ H(x,Dlo,uj) > 

(35) 

where H(-) = is the HJB equation for the 
problem. 



Note that these assumptions hold for the particular 
choice of discretization via transition probabilities 
that we have outlined. Under these assumptions we 
proceed to look at convergence results for certain 
approximations to the value function. 

5.2. Convergence of the Approximation 

In this section we prove the validity of the ap- 
proximations to the HJB. Our aim is to prove the 
convergence of the two terms defined below to the 
viscosity solution. 

V*{x):=Yimsn V V h {y), (36) 
Ho 

V*(x):= lim inf ^(y). (37) 

y 

HO 

The following is a generalization of [23, Ch.IX, 
Lemma 4.1, Theorem 4.1] 

Lemma 14. V* is a viscosity sub-solution to the 
HJB H(x,DV,V) = 0. 

Proof. Let ui e C^L) s.t V* - u> has a local 
maximum at x € L. Note that without loss of gen- 
erality we can assume that x is a strict maximum 
on L by redefining uj suitably [22]. There exists 
some subsequence h converging to zero s.t V h — uj 
has a maximum at yh on L^ (such that yh — > x as 
h goes to zero). From the local maximum property 
it follows that 

V h (y h ) - uj(y h ) > V h (y) - u>(y) , y G L° h . (38) 

Hence taking a 7^ € (0, 1] we have 

7,; 1 (V h (y h ) - u(y h )) > V h (y) - oj(y) , y e K 

(39) 

=> u(y) - 7fc _1 w(|/ fc ) > V h (y) - 7fc _1 V h (y h ). 

(40) 

Applying the operator F h and using the properties 
that it satisfies we obtain 

F h (oj)(y h ) - uj(y h ) > F h (V h )(y h ) - V h (y h ). (41) 

At y = yh the right hand side of the above equation 
is zero (from Eq(34)). Dividing by h and taking the 
limit as h tends to zero, we have 

lim ^M^>0. 
h 

hi.0 

,. F h (oj)(x) -u(x) n 

lim v A ; ^ > (42 

h-tO h 
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Using Eq (35) we obtain 

H(x,Duj 1 uj) < 0. 
Hence the theorem is proved. 

Similarly we can show that V* is a viscosity super- 
solution to the HJB equation. We then have the 
following result. 

Theorem 2. Assume that V h converges to the 
boundary conditions in a uniform manner i.e 

limV h (y)=V(x), xedT. (43) 
Ho 

We then have 

lim V h (y) = V(x), x e L. 

hlO 

Proof. V* was shown to be a viscosity super- 
solution to the HJB. Hence using Eq (43) we have 
from standard comparison theorems for continuous 
viscosity solutions [17] that 

V(x) <V,(x), x e L. (44) 

Similarly since it was shown that V* is a viscosity 
sub-solution to the HJB, we apply the comparison 
theorem to obtain 

V(x) >V*(x), x e L. (45) 

From Eqns (44), (45) we have that 

V* < V* on L. 

However from the definitions of V* and V* we have 
that 

V* > V, on L. 

Hence using the two inequalities above we have 
V* = 14 on L and the statement of the theorem 
follows. 

6. Simulations 

We now provide an example problem on SU (2) to 
demonstrate the results of this article. The system 
dynamics is given by 

U =[v 1 I x +v 2 I z ]U, \\v\\=2 

(46) 

U(0) = U , 
U ,U e SU(2). 



The value function for this problem is given by 

tu (v) 

S{U ) = inf / exp{-A s}ds. (47) 



As previously mentioned, cxp(g) generates G. 
Hence instead of obtaining several complexes, 
patching them together and refining them on the 
areas of overlap, we take a sufficiently large area 
of the algebra su(2) and perform a triangulation 
on it. This is followed by 'glueing' together the 
value function at points corresponding to the iden- 
tity element of SU(2). This simplified discretiza- 
tion will yield accurate solutions close to the target 
set (and on the interior of the region being con- 
sidered), but loses accuracy towards the end of the 
region. The mapping exp(x I x + y I y + z I z ) from 
su(2) to SU(2) provides three natural parameters 
with which to visualize the value function on the 
Lie algebra. Note that the map is not injective and 
multiple points from the algebra may be mapped 
to the same point on the group. 



0. 




Figure 3: Normalized minimum time function with a control 
of norm 2. 

We use standard numerical tools from compu- 
tational geometry to obtain a triangulation of the 
chosen subset of the algebra visualized as a region 
in M 3 (since su(2) is isomorphic to M 3 ). Perform- 
ing the value iteration mentioned in Section 5 the 
resulting value function is shown in Fig 3. Lighter 
shading at a point indicates a larger value of the 
minimum time function at that location. The stop- 
ping criteria used in the algorithm designed is as 
follows. At each step of the iteration the abso- 
lute change in the value function at all points on 
the mesh is computed. The maximum value of this 
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Figure 4: Optimal trajectories to the identity element start- 
ing from different points on each of the axes. 



E o.o 
1 




Iteration count 



Figure 5: Maximum of the absolute change in the value func- 
tion over the grid at each iteration. 



change across the entire mesh is determined and a 
threshold value for this stopping metric, at which 
the algorithm should terminate, is set. This metric 
is indicated in Figure 5. Note that there are several 
different possibilities choices for a stopping crite- 
ria. We defer to future work the analysis of various 
possible stopping metrics and a detailed quantita- 
tive study of numerical convergence rates and error 
bounds. 

The optimal controls are obtained using a dis- 
cretized version of the verification theorem and 
from them we can obtain sample trajectories from 
various starting points as indicated in Fig 4. Note 



that this figure must be carefully interpreted since 
the mapping from this region into ST/ (2) is not in- 
ject ive. 



7. Conclusion 



In this paper we introduced a rigorous frame- 
work for the numerical techniques involved in using 
the Dynamic Programming technique from optimal 
control theory for the control of quantum spin sys- 
tems evolving on compact Lie groups. Numerical 
simulations were preformed by triangulation of the 
group which, due to the well studied numerical pro- 
cedures available for tesselation of surfaces, enable 
better numerical speed and efficiency of implemen- 
tation. In addition, the solution can be made more 
accurate at points of the group where such accuracy 
is desired (e.g around the origin in Fig. 3 where the 
solution is non differentiable). The dynamic pro- 
gramming methods provide a framework that can, 
in principle, be used for systems with an arbitrary 
number of qubits unlike limitations on the Lie the- 
oretic methods. In addition, alternative numerical 
techniques that use the calculus of variations are 
subjected to issues in the entrapment at local min- 
ima - a drawback absent in the current approach. 

The value function iteration methods when used 
on any grid, suffer from the curse of dimensional- 
ity and hence become intractable for higher dimen- 
sional systems. For instance, the number of spatial 
dimensions in a quantum spin 1/2 system with n 
qubits grows as 4" — 1. Possible directions of future 
work may involve a study of methods such as fast 
marching [24] or meshless techniques to improve 
the speed of computations. Inspired by the dy- 
namic programming framework in this article and 
the curse of dimensionality free approaches in [25] , 
new methods are currently being developed for re- 
duced dimensionality approximation techniques to 
quantum control. 

There exist classes of control problems such as 
those involving quantum systems with bounded 
controls and drift for which the value function is dis- 
continuous. Viscosity solution techniques for such 
discontinuous cost functions may be used to pro- 
vide the technical framework for the use of impul- 
sive controls in the dynamic programming approach 
to quantum control. 
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